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Abstract. We consider a generic type of dark energy fluid, characterised by a 
constant equation of state parameter w and sound speed c s , and investigate the 
impact of dark energy clustering on cosmic structure formation using the spherical 
collapse model. Along the way, we also discuss in detail the evolution of dark energy 
perturbations in the linear regime. We find that the introduction of a finite sound speed 
into the picture necessarily induces a scale-dependence in the dark energy clustering, 
which in turn affects the dynamics of the spherical collapse in a scale-dependent way. 
As with other, more conventional fluids, we can define a Jeans scale for the dark energy 
clustering, and hence a Jeans mass Mj for the dark matter which feels the effect of 
dark energy clustering via gravitational interactions. For bound objects (halos) with 
masses M 3> Mj, the effect of dark energy clustering is maximal. For those with 
M <C Mj, the dark energy component is effectively homogeneous, and its role in the 
formation of these structures is reduced to its effects on the Hubble expansion rate. 
To compute quantitatively the virial density and the linearly extrapolated threshold 
density, we use a quasi-linear approach which is expected to be valid up to around the 
Jeans mass. We find an interesting dependence of these quantities on the halo mass 
M, given some w and c s . The dependence is the strongest for masses lying in the 
vicinity of M ~ Mj. Observing this M-dependence will be a tell-tale sign that dark 
energy is dynamic, and a great leap towards pinning down its clustering properties. 
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1. Introduction 

The apparent accelerating expansion of our universe is, according to the standard model 
of cosmology, best described by the presence of a dark energy component with a strong 
negative pressure, dominating the gravitational physics on large scales (see, e.g., [1-3] 
for reviews). However, our knowledge of the actual properties of this dark energy is very 
limited, and a number of open questions remain — Is dark energy dynamic or not? Does 
it take part in clustering? Is the dark energy's behaviour scale-dependent? — to name 
but a few. 

In order to distinguish between various models of dark energy, we must appeal to 
a variety of observational tests, each probing a different aspect of the dark energy's 
dynamics. Luminosity and angular diameter distance measurements using, respectively, 
type la supernovae and the baryon acoustic oscillation scale have yielded — and will 
continue to yield — interesting information on the dark energy's influence on the 
expansion history of the universe. The cosmic microwave background temperature 
anisotropies [4-8], as well as their cross-correlations with tracers of the large-scale 
structure distribution [9-12], provide a means to track the dark energy's impact on 
the evolution of the gravitational potential via the late integrated Sachs- Wolfe effect. 
Weak gravitational lensing of distant objects probes the dark energy's effect on the 
distance-redshift relation and the growth function [13-15]. Lastly, the formation of 
gravitationally bound objects such as galaxies and galaxy clusters is also sensitive to 
the detailed properties of the dark energy component [16-19]. 

Into the last category falls the so-called spherical collapse model [20], which, as 
the name suggests, is a model of gravitational collapse simplified by the assumption 
of spherical symmetry. In the model, a spherically symmetric overdense region with 
uniform density evolves to a configuration of infinite density under its own gravity, and 
a gravitationally bound object is said to be formed. 

The original spherical collapse model was constructed under the assumption of 
an Einstein-de Sitter (EdS) universe. Already by the 1980s and the early 1990s the 
model had been extended to include a cosmological constant [21-23], and later for 
quintessence [24]. More analyses have followed since then, and all reached the conclusion 
that dark energy has an important impact on the formation of gravitationally bound 
structures [25-28] . The topic appears to have gained momentum again during the past 
year [29-34] mainly because of the positive expectation that these results will be testable 
against observations in the not-too-distant future [35]. 

In this paper we investigate by means of the spherical collapse model how a 
generic dark energy component characterised by a constant equation of state parameter 
w < —1/3 and sound speed c s affects the formation of cosmic structures. In previous 
studies the dark energy component is generally allowed to take on different values of w. 
However, the sound speed c s is invariably assumed to be either approaching the speed 
of light so that the dark energy is essentially homogeneous, or exactly vanishing so that 
the dark energy fluid is comoving with the matter component. These assumptions 
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undoubtedly simplify the calculations considerably and represent the two limiting 
behaviours of dark energy clustering. However, the case of a general sound speed remains 
interesting in that it introduces a scale-dependence into the problem. Identifying this 
dependence will give us yet another pointer to the true nature of dark energy. Here, 
we demonstrate for the first time how to incorporate a dark energy component with an 
arbitrary sound speed into the spherical collapse model (in an approximate way). 

The paper is organised as follows. In section 2 we introduce the evolution equations 
for the spherical collapse and the corresponding equations of motion for the dark energy 
component. In section 3 we discuss dark energy clustering within the framework of 
linear perturbation theory. Section 4 contains our numerical results. Our conclusions 
are presented in section 5. 



2. Spherical collapse 

2.1. The spherical top hat and equations of motion for the matter component 

In its most basic formulation, the spherical collapse model assumes there exists a 
spherically symmetric overdense region on top of an otherwise uniform background 
matter density field. The overdense region is characterised by a physical radius 
Ri = R{ti) at the initial (conformal) time T;, and a uniform initial energy density 

where p m (T) denotes the energy density of the background matter field. This is our 

spherical "top hat" perturbation, and the mass contained within is given by 

Ait Att 
M = -p m {n){l + CM = yPm(ro)(l + C,)Xl (2.2) 

where tq denotes the present time, and we have defined 

X=- (2.3) 
a 

as the comoving radius of the top hat. 

The evolution of the physical top hat radius R with respect to cosmic time t is 
described by the familiar equation of motion 

where we have incorporated in the equation the possibility of a second energy component 
with a nonzero pressure denoted by the subscript Q. This second component is uniform 
inside the top hat region defined by the radius R, and shall be our dark energy 
component in this work. Equation (2.4) can be equivalently expressed as an equation 
of motion for for the comoving top hat radius X with respect to conformal time r, 

§ + ?4 = -^r*M + no- + ^Ki (2.5) 



Spherical collapse of dark energy with an arbitrary sound speed 



4 



where • 



d/dr, H = aH is the conformal Hubble parameter, and 

o SPr 



(2.6) 



rest 



is the square of dark energy sound speed defined in the dark energy fluid's rest frame. 
Note that in identifying c s in equations (2.5) and (2.6) as the rest frame sound speed, 
we have implicitly assumed that we are dealing only with length scales much smaller 
than the Hubble length, and that averaged over the spherical region, there is no bulk 
flow of dark energy relative to the dark matter fluid. We also assume c 2 s to be constant 
in time and space. 

Since the total mass of matter inside the top hat M = (4n /3)p^R 3 is conserved, 
the top hat matter density contrast 5^ can be easily expressed as a function of the top 
hat radius, 



S. 



tlw 



Pm{r) 



- 1 



'1 + 6. 



th 



o(t) Rj 
R(t) 



a(Ti 



-l = (l + « 



th -i 
m,iJ 



X{T) 



- 1. 



(2.7) 



For the dark energy density contrast Sq , two limiting cases have been studied in the 
literature. The first is the non-clustering limit, in which the dark energy sound speed c s 
is taken to approach the speed of light, see e.g. [24]. In this case, Sq is effectively zero, 
so that the role of dark energy in the spherical collapse enters only through the Hubble 
expansion of the background.! 

The second is the "co moving" or clustering limit, in which the dark energy sound 
speed is exactly zero, see e.g. [31]. As we shall see in the next section, the Euler 
equation for the dark energy fluid in this limit is identical to its counterpart for a 
nonrelativistic dark matter fluid. This means that the bulk velocity fields of the two 
fluids are the same; the fluids are thus said to be comoving. Note that this observation 
does not imply the dark energy and dark matter density contrasts evolve in the same 
manner, since the conditions for energy conservation differ between the two fluids. It 
does, however, imply a conservation law for the dark energy component inside the top 
hat, so that the evolution of pg can be simply expressed as 



dp% 3 dR 



dt 



Pc 



0. 



or 




+ Pc 







(2- 



(2.9) 



in terms of comoving quantities. 

Strictly speaking, these two limiting cases are the only ones for which the top 
hat formulation is exact. The case of an arbitrary dark energy sound speed c s is 

| Note that in order to get exactly zero dark energy clustering, the sound speed would have to be 
infinite; See the linear solutions (3.4) and (3.6). 
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strictly not amenable to this simple treatment, since the existence of a finite sound 
speed and therefore the provision for the propagation of sound waves imply that the 
energy densities — both dark matter and dark energy — inside the overdense region must 
evolve to a nonuniform configuration, even if they are initially uniform. Having said this, 
however, we must also bear in mind that the spherical collapse model is itself a simplified 
model of structure formation, and the top hat density contrast should be interpreted as 
the average density contrast inside a region after a unit top hat filtering function has 
been applied. If we take this as our guiding principle, then the generalisation of the 
spherical top hat to include a dark energy component with an arbitrary sound speed 
simply requires that we interpret <5q as the spatially averaged density contrast of the 
dark energy field inside a region of comoving radius X. Symbolically, this spatial average 
can be expressed as 

^ h (r) = ^ j* dx xH q {x,t), (2.10) 

where x = \x\, and x denotes the comoving coordinates. 



2.2. Equations of motion for the dark energy component 

It remains to specify an evolution equation for the dark energy density perturbation 
Sq(x, t). We begin by writing down the continuity and Euler equations for a relativistic 
fluid a in an expanding background in the pseudo- Newtonian approach [36], 

pa + m{p a + P a ) + V ■ [(p a + P a )u a ] = 0, 

VP + u P 

ii a + Uu a + (u a -V)u a + ^ + V^ = 0. (2.11) 

Pa + Pa 

Here, V = d/dx, u a is the peculiar velocity of the fluid, and the potential 4>n can be 
obtained from the Poisson equation 

V 2 <j) = 47rGa 2 J2$Pa + 35P a . (2.12) 

a 

These equations should apply if we restrict our considerations to (i) length scales 
much smaller than the Hubble length, (ii) nonrelativistic peculiar velocities, and (iii) 
nonrelativistic sound speeds c, < 1. We demonstrate in Appendix A that, at the 
linear level, these equations are indeed consistent with the Newtonian limit of a general 
relativistic formulation (see, e.g., [37-39]). 

Defining the equation of state parameter for the dark energy component 

w = ^-, (2.13) 
PQ 

and using the definitions P a = P a + SP a and 5^ = p~ 1 5P Q ,, equation (2.11) can be 
rewritten for a = Q as 

5 Q + 3H(6% - wS Q ) + V • \{p Q + P Q )u Q /p Q ] = 0, 

u Q + Uu Q + (u Q ■ V)u Q + — ^— v + V<p N = 0. 2.14 

1 + w + dn + bh 
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Since in our set-up the universe contains only dark matter and dark energy, the Poisson 
equation (2.12) now reads 

V 2 = A7rGa 2 [p m 5 m + p Q (5 Q + 35 J)], (2.15) 

where, for our particular problem, the dark matter density perturbation 5 m (x, r) takes 
the form 

6 (xr) = { ^ (r) ' (2 16) 

° m[Xl ) \ 0, x>X(r). [ ' 

with the top hat density contrast 5J£(t) given by equation (2.7). 

The continuity and Euler equations (2.14) are nonlinear in the quantities Sq and 
uq, which even under the assumption of spherical symmetry are nontrivial to solve. 
Therefore, as a first approximation, we linearise them to obtain 

^ n + 3H(5j lin - w5 l £) + (1 + w)9l n = 0, 

V 2 <5 p,lin 

6 l £ + (1 - 3w)H6 l X n + - Q + V 2 <f) N = 0, (2.17) 

1 + w 

where we have defined the divergence of the dark energy velocity field to be 9q = V-Mq, 
and assumed w to be constant in time. Linearisation assumes that the perturbed 
quantities Sq and 9q are small. This is likely a good assumption since (i) the presence of 
a finite sound speed c s naturally hinders the clustering of dark energy, keeping Sq ~ Sq 1 
small relative to S m , and (ii) even in the limit c s = where the clustering of dark energy is 
most efficient, a fully linear analysis shows that Sq 1 is suppressed relative to because 
of the dark energy's negative equation of state parameter (see section 3). Either way, 
the assumption of linearity in dark energy clustering can be easily checked a posteriori 
against solutions of the evolution equations for consistency. Finally, let us stress again 
that we are linearising only the dark energy equations of motion; the evolution of the 
dark matter component is still fully nonlinear, and described by the spherical collapse 
detailed in the previous section. We shall call this the "quasi-nonlinear" approach. 

Upon linearisation, it is useful to recast the equations of motion in Fourier space. 
Define the Fourier transform for some field A(x,t) as 



A(x, t) = - — / d 3 k A(k, t) exp(ik ■ x) 
(27rr J 



(2k) 

1 f dk k 2 A(k,r)^^. (2.18) 



2ir 2 J ' kx 

Then, using the Poisson equation (2.15) and the relation <5y>' lm = c^Sq 1 + 37^(1 + w)(c 2 s — 
w)0Q 1 /k 2 [6], equation (2.17) can be equivalently expressed as 

5 Q + 3(c 2 - w)H5^ + (1 + w)9% = 0, (2.19) 

Jin „ k 2 C 2 ~ ~ ~ 

9 Q +(1- 3c 2 s )U9% - i^hJ? + ^Ga 2 [p m 5 m + p Q (l + 3c 2 s )8%] = 0, 

where we have dropped subdominant terms proportional to / H 2 /k 2 , since we are 
interested only in subhorizon scales k ^> %. For the dark matter density contrast 
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S m (x,r) given in equation (2.16), the Fourier space equivalent is 



5m (fc, t) = Att dx x 2 5%(t) 
Jo 



kx 



Att Att 
5t{r)-X*W{kX) = y[(l + - X 3 ]W(kX), (2.20) 



where 

(kX) 



W(kX) = Ji: ^[sm(kX) - kX cos(kX)} (2.21) 



by convention. 

Lastly, we would like to relate b^qik, r) to the average dark energy density contrast 
inside the top hat, Sq-(t), as defined in equation (2.10), since this is the quantity that 
ultimately governs the evolution of the top hat radius X via equation (2.5). This step is 
simple: we only need to identify Sq(x, t) with 5g n (x, r), the latter of which is obtained 
by Fourier transforming 5g n (/i;,r). Thus, equation (2.10) simplifies to 

6 Q M = ^J dk k 2 W(kX)~5 l ^(k, t). (2.22) 
Our set of equations of motion is now complete. 



3. Linear theory 

Before we present the results of the spherical collapse model, let us first consider the 
evolution of dark matter and dark energy perturbations in the linear regime, i.e., where 
the dark matter perturbations are also tracked with linearised equations of motion. This 
exercise is useful for two reasons. Firstly, as we shall see, an understanding of the linear 
evolution can shed light on many essential features of the dependence of dark energy 
clustering on its equation of state parameter w and sound speed c s . At the same time, 
the linear solution also sets the initial conditions for the spherical collapse model. 

Secondly, some semi-analytic theories of structure formation such as the Press- 
Schechter formalism [40] and the excursion set theory [41-44] require as an input a 
linear critical density contrast 5'™ n . In these theories a collapsed structure is assumed 
to have formed once the linearly evolved matter density contrast reaches the threshold 
value 5 l ™ n at some time r co n. The value of t^Si can De determined from the spherical 
collapse model by interpreting r co n as the instant at which the top hat radius vanishes. 
In practice, this means that in order to extract 8^ for a particular cosmological model, 
we need to solve both the nonlinear and the linear equations of motion at the same time. 

We have already written down the linearised equations of motion for the dark 
energy perturbations in equations (2.17) and (2.19). For the dark matter component, 
the corresponding equations are 

e + + 47rGa 2 [p m C + PqO- + 3c ^q] = °- t 3 - 1 ) 
These and equation (2.19) are solved simultaneously to determine the evolution of 5^ 
and #l n . 
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3.1. Dark energy evolution in the linear regime 

Let us consider first the linear evolution of the dark energy perturbations. Here, it is 
convenient to combine the first order differential equations (2.19) for Sq 1 and 8q n into one 
second order differential equation for Sq 1 , and also adopt a new time variable s = In a. 
Assuming w and c s to be constant in time, we find 

S Q n " + V(s)~6 Q »' ■ — - 3 

where ' = d/ds, and 

V(s 



1,2 2 



q =^l+w)n m (s)6^, (3.2) 




«(s) = 3w [ 1 + — J + -(1 + w)n Q ( 8 ), 
X(s) = 1 + 3^ 1 + ^ - 3(c 2 - «,) - 2(1 + «,)n (a) 



(3.3) 



For the cosmological models considered in this work, T>(s) and \k(s)\ are of order unity, 
while X(s) ~ 1 always holds true because of our assumption of k ^> "H. 

Equation (3.2) describes a damped harmonic oscillator with a driving force sourced 
by perturbations in the dark matter fluid. Exact analytic solutions do not exist for 
arbitrary cosmologies. However, approximate solutions can be constructed in certain 
limits: 

(i) Clustering limit. This is the limit in which k 2 c1/'H 2 <C |k| ~ 1. In this case, all 
coefficients in the differential equation are of order unity. It is therefore necessary 
to specify the exact time dependence of T>(s), k(s) as well as in order to find 
a solution. Formally setting k = 0, the solution is particularly simple during the 
matter domination epoch, where fi m (s) ~ 1, J)q(s) 1, %' ' /% — —1/2, oc a 
and 6*^ n ~ — H5]" 1 . At s — s» ^> 1 it has the asymptotic form 

?lin _ 1 + w Bin 

d Q " 1-3^^' 

-^[3(^-^ + 1]^, (3.4) 

where we have obtained the solution for 6q n from the continuity equation by first 
differentiating Sq 1 with respect to time. 

At first glance, the solution (3.4) for Sq 1 appears to be at odds with the solution 
obtained in, e.g., reference [39] in the same limit (i.e., < I; < K s , where 
% s = H/c s is the inverse of the sound horizon, or the "Jeans wavenumber" kj as we 
define in equation (3.5) below). In particular, the solution of [39] depends explicitly 
on the sound speed c 2 , whereas our solution does not. Part of the discrepancy can 
be traced to the term X{s) defined in equation (3.3). In our analysis we always 
approximate this term as X{s) = 1, while some contributions proportional to "H 2 //c 2 
have been retained in the analysis of [39]. 



Spherical collapse of dark energy with an arbitrary sound speed 



9 




However, we believe that this discrepancy is of little consequence. As we 
demonstrate in Appendix B, the H <C k <C % s limit is well-defined only for those 
dark energy sound speeds satisfying c 2 s < 10 -3 . Thus, from a numerical point of 
view, our approximate solution and that of [39], where they are actually applicable, 
are consistent with one another to better than one part in a thousand. 

(ii) Non- clustering limit. This limit corresponds to k 2 c 2 /1-L 2 ^> |ac| ~ 1, which is also 
the steady-state limit (\8q" /8q\, \"q I& X q\ 1). The solution can be obtained by 
formally setting 5 X q" = Oq' = k = 0. Defining the "Jeans wavenumber" 

kj = -, (3.5) 

C s 

the steady-state/non-clustering solution then reads 

3 
2 

^ -f^[cl-wn m {s)]5^. (3.6) 

Note that, unlike the clustering solution (3.4), the non-clustering solution is not 
restricted to the matter domination epoch. Furthermore, the derivation of (3.6) 
does not in fact require the assumption of a linear 8 m , since 8 m enters into the 
differential equation (3.2) only through the gravitational potential 0, and hence the 
Poisson equation, which is in any case linear in 8 m . This means that the steady- 
state/non-clustering solution (3.6) would have been equally valid had we replaced 
5g n with the Fourier transform of the top hat density contrast 5^. We shall make 
use of this solution again later on in the analysis. 

The form of the non-clustering solution is akin to those commonly found in hot or 
warm dark matter scenarios, in which kj is associated with the free-streaming scale 
of the problem (see, e.g., [45]). However, since dark energy has a non- zero w while 
free-streaming dark matter does not, an extra pref actor (1 + w) is incurred in the 
solution (3.6). 

(iii) Unstable limit. So far we have implicitly assumed c s to be a real number. Let 
us entertain ourselves for a moment with the possibility of an imaginary dark 
energy sound speed. In the limit \k 2 c 2 /'H 2 \ <C 1, the dark energy perturbations 
are described by the same clustering solution as equation (3.4). Contrastingly, the 
\k 2 c 2 /'H 2 \ ^> 1 limit is unstable. Formally setting T> — k — = 0, equation (3.2) is 
solved in the matter domination epoch by Sq 1 = CiIo(uy/a) + ^^(wv^li where 
Iq(x) and K (x) are the zeroth order modified Bessel functions of the first and 



the second kind respectively, and u = 2k\c s \/y H$VL m . For x = u^fa ^> 1, the 
modified Bessel functions have the asymptotic forms Iq(x) ~ exp(x) /y/2Tix and 
Kq(x) ~ ^7i/2xexp(—x). Thus, we find for the linear dark energy density contrast 
the asymptotic solution 

^ n ~a- 1 / 4 exp(o;v / a). (3.7) 
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This exponential growth of 5 1 q in turn sources the evolution of the dark matter 
density contrast via the Poisson equation (2.15). Consequently, 5 l ™ also exhibits a 
similarly explosive and strongly scale-dependent growth that at first glance appears 
to be in conflict with our current understanding of large-scale structure formation 
unless \c s \ is very small. We shall therefore not pursue the case of an imaginary 
dark energy sound speed any further in the present work. 

Given the limiting solutions (3.4) and (3.6), we can try to interpolate between the 
clustering and non-clustering regimes using the following (rough) interpolation formulae: 

?lin _ 1 + w ~ lin 

~ l-3w + (2/3){k/kj) 2 m ' 



aim ^ 

Q = 



1 + W 



of 2 \ 1-3W 

3{c s - w) 



(3-8) 



l-3w + (2/3)(k/kj)\ 

The maximum error is 30% at k ~ kj. These interpolation formulae are valid during 
the matter domination regime, and can be used to set the initial conditions for the dark 
energy component in the spherical collapse model. § 

Lastly, let us define a "Jeans mass" scale analogous to the Jeans wavenumber kj 
given in (3.5), i.e., the mass scale at which we expect the effects of the dark energy 
sound speed to set in. The Jeans mass is defined here as 

, . Air . . / Xj(a) \ 3 
Mj(a) = — p m (a) , (3.9) 



3 ' y ' \ 2 

where Aj ee 2n/kj. Evaluating the expression at a = 1, we find 

Mj = 9.7 x 10 23 n m c 3 H^Mq. (3.10) 

For example, given Q m = 0.3 and h = 0.7, we have Mj = 1.3 x 1O 16 M for c 2 = 10~ 5 , 
and Mj = 4 x 1O 14 M for c 2 = 10 -6 today. Note that the mass here refers to the mass 
of the dark matter component, not the dark energy! 



3.2. Linear threshold density 

Our second motivation for considering linear theory is the computation of the linear 
threshold density, defined as 

<di = C' lm ( r coii), (3.11) 

where S^ ,hn (r) is the linearly evolved top hat matter density, and r co u is the instant 
at which the top hat radius goes to zero. As the name implies, S^' hn (r) is the linear 
version of the quantity S^(r) defined in equation (2.16), and is tracked by the equations 

§ The full equation (3.2) in fact has an exact analytic solution encompassing all three limits discussed 
above in terms of Bcsscl functions in the matter domination epoch. However, the complexity of the 
solution rather obscures the simple physics behind the problem. We therefore do not quote it here. 
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of motion (3.1) upon the replacements 

jlin/™ ^\ i /ith,lin 



Sq(x,t) -+ <^ lin (r), (3.12) 



where 



^ lin (r) = ^- 2 Jdk k 2 W(kX)~5^(k,r) (3.13) 
is the linearly evolved dark energy density contrast averaged over the top hat volume. 

4. Numerical results 

In this section we proceed to solve numerically the evolution equations for the spherical 
collapse model presented earlier in section 2.|| We assume a flat spatial geometry for the 
universe so that the dark energy fraction today is related to the dark matter fraction by 
Qq = 1 — Q m . For the choice of the parameter fl m and the present Hubble rate H , we 
use the WMAP 7-year best-fit values [46]. We consider only those cases with constant 
w and c s , although our formulation is applicable also to scenarios with time- dependent 
dark energy parameters. 

We begin the evolution at a dimensionless time coordinate of 

UH = 2 x 10" 6 , (4.1) 

corresponding to an initial scale factor of 

2/3 



<h = ao I ^ j ' ( 4 ' 2 ) 

if we assume t{ to lie well within the matter domination epoch. Taking an initial matter 
overdensity of 8 t ^ li the initial value of the top hat radius can then be obtained directly 
from equation (2.2) given some mass M. This mass, which we dub the "halo mass", is 
also the mass of dark matter contained in the final collapsed object. % Unless otherwise 
stated, the initial matter density contrast is taken to be 

5% = 3 x 10- 4 . (4.3) 

We have chosen the above values for the initial time and matter density contrast so that 
the collapse occurs at a time when the dark energy component constitutes a significant 
part of the universe's energy budget. 

Since the top hat evolution equation is a second order differential equation, we 
must also specify the time derivative of R. This can be constructed by differentiating 

|| The numerical code is written in C++ and employs GNU scientific libraries for solving the evolution 
equations and for interpolating the integral in equation (2.22). A verification of the convergence of 
equation (2.22) against the number and spacing of /c-bins has been conducted. 

% Note that the mass M here refers to the mass of the dark matter component only, although as 
suggested in [31], the true mass of the bound object should in principle include the contribution from 
the clustered dark energy component as well. 



Spherical collapse of dark energy with an arbitrary sound speed 



12 



equation (2.7) with respect to time. Because the initial matter density contrast is 
much less than unity, we can approximate d5^/dt ~ d5^' hn /dt ~ H5^ ,lm using linear 
perturbation theory, and thus, 



Finally the initial conditions for the dark energy evolution is given in section 3.1, 
particularly by the interpolation formula (3.8). 

4-1. The collapse 

Figures 1 to 4 shows the physical top hat radius normalised to the initial radius as 
a function of the dimensionless time coordinate tHo for several choices of w, c s and 
halo masses M. The corresponding matter overdensity (2.7) and the dark energy 
overdensity (2.22) are also shown in juxtaposition. 

In figure 1 we present the results for a dark energy component with c 2 s = 10 _1 and a 
halo mass of M = 10 14 M Q for various equation of state parameters. These choices of c s 
and M satisfy M <C Mj according to equation (3.10), and ensure that we are in the non- 
clustering regime. Since dark energy clustering is minimal (< 10~ 5 relative to the dark 
matter density contrast), only the equation of state parameter w, i.e., the homogeneous 
part of the dark energy fluid, plays a role in the dynamics of the spherical collapse. 
Indeed, as we see in the figure, the less negative the equation of state parameter, the 
later the collapse. The reason is that for the same Qq, the less negative w is, the earlier 
the dark energy comes to dominate the energy content of the universe, thereby inhibiting 
the growth of the top hat overdensity through Hubble expansion from an earlier stage. 
Comparing the ACDM case and a model with w = —0.7, the collapse time is delayed 
by some 20%. For equation of state parameters more negative than —1, the opposite 
trend is seen; in the case of w = — 1.3 the collapse occurs some 10% faster than in the 
ACDM limit. 

In figure 2 we fix M = 1O 16 M , but vary the dark energy sound speed and equation 
of state. We choose w = —0.7 and w = —0.8, since, based on results from linear theory 
(see section 3), dark energy clustering is most enhanced by a deviation of w from — 1 
in the positive direction. Although the effect of dark energy clustering on the spherical 
collapse is quite small, the trend is clear: the smaller the sound speed, the faster the 
collapse. This is to be expected, since the smaller the sound speed is, the more efficiently 
the dark energy component clusters, and this clustering in turn contributes to sourcing 
the collapse of the top hat on the r.h.s. of equation (2.5). Note that although the dark 
energy component exhibits some degree of clustering in these cases, the density contrast 
for almost the entire collapse history is quite small until the last moments when R — > 0. 
This indicates that our quasi-nonlinear approach — in which the dark energy component 
is evolved with linearised equations — is valid for the model parameters adopted in this 
figure. 




(4.4) 
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Figure 1. Top: Spherical top hat radius R normalised to the initial radius Ri 
as a function of the dimensionlcss time coordinate tHo. The reference EdS and 
ACDM models are represented by the blue dot-dash and the red dot-dot-dot-dash 
lines respectively. The two black solid lines denote, from right to left, the cases of 
w = —0.7, —0.9, while the two green dashed lines denote, from right to left, the cases 
oi w = — 1.1,— 1.3. For these four cases, we have chosen the dark energy sound speed 
to be c 2 s = 10 _1 , and a halo mass of M = 10 14 M Q . Bottom: The corresponding top hat 
matter and dark energy density contrasts. For the cases of w = —1.1, —1.3, the dark 
energy density contrasts are negative, i.e., they are underdensities, and the 5q values 
presented in this plot are absolute values. Note that in the reference EdS and ACDM 
models, there is no dark energy clustering. For all cases the initial matter overdensity 
has been chosen to be i5^ i = 3 x 10~ 4 . 

Figure 3 shows the cases of a fixed sound speed c 2 s = 10~ 4 , equation of state 
parameters w = —0.7, —0.8, and three different halo masses. Figure 4 is similar, but 
with the sound speed fixed at c 2 s = 10~ 6 . For these sound speeds, the corresponding 
Jeans masses Mj are 4 x 1O 17 M and 4 x 1O 14 M respectively. In both figures, we see 
that the larger the halo mass, the faster the collapse. This can be understood from the 
non-clustering solution (3.6) (valid here since (almost) all halo masses considered are 
less than the Jeans mass). Since the dominant Fourier mode is that corresponding to the 
comoving top hat radius X which is itself associated with the halo mass M, a reasonable 
generalisation of the non-clustering solution (3.6) for the dark energy component in 
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Figure 2. Same as figure 1, but for a fixed halo mass M = 10 16 Mq. The black 
solid lines indicate w = —0.7 and, from right to left, a dark energy sound speed of 
c 2 s = 10 -1 ,10 ,10 -6 . The green dashed lines denote w = —0.8 for the same set of 
dark energy sound speeds. The reference EdS and ACDM models are again represented 
by, respectively, the blue dot-dash and the red dot-dot-dot-dash lines. The right panels 
show the same results, but zoomed in on the time interval tHo = 1.05 — > 1.25. 

terms of the halo mass would be 




The expression clearly shows that for a given sound speed, the absolute value of the 
dark energy density contrast increases with halo mass M. The enhanced dark energy 
density contrast in turn hastens the collapse of the dark matter top hat. 

Note that for the case of M = 1O 16 M and c\ = 10 -6 , the dark energy density 
contrast is of order 0.1 for much of the collapse history. This suggests that our 
approximation scheme for the dark energy evolution is approaching its limits of validity; 
the approximation breaks down for larger masses. Thus the rule of thumb regarding the 
quasi-nonlinear approach appears to be that it can be safely used for halo masses up to 
roughly the Jeans mass Mj, but not beyond. + 

+ Obviously, we came to this conclusion based on the rather extreme cases of w = —0.7 and w = —0.8. 
In general, however, we expect the validity of the quasi- nonlinear approximation scheme to be dependent 
also on the choice of w, where the less w deviates from —1, the higher the halo mass for which the 
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Figure 3. Same as figure 2, but for a fixed sound speed cj. = 10~ 4 . The black solid 
lines denote w = -0.7 and, from right to left, halo masses of M = 10 12 , 10 14 , 10 16 M o . 
The green dashed lines denote w = —0.8 for the same set of halo masses. The blue dot- 
dash and the red dot-dot-dot-dash lines represent the reference EdS and ACDM models 
respectively. The right panels show the same results in the time interval between 
tH = 1.05 and tH = 1.25. 



4-2. Linear threshold density 

Next, we compute the linear critical density contrast <5j.™ u required for use with such semi- 
analytic theories as the Press-Schechter formalism and the excursion set theory. This 
can be achieved by solving simultaneously both the nonlinear and the linear equations 
of motion for the spherical collapse, and formally identifying 8^ as the linearly evolved 
matter density contrast at the instant the top hat radius vanishes. Figure 5 shows 5 l ™ a 
as a function of the halo mass for various combinations of w and c s . The initial matter 
overdensities are chosen such that all halos collapse at z = (top panel) , z — 1 (middle) , 
and z = 2 (bottom). 

Clearly, for the reference cases of an EdS and a ACDM universe, the linear threshold 
density is independent of the halo mass. However, once a finite dark energy sound 
speed is introduced into the picture, 5 l ™ n becomes mass- dependent, with 5 l ™ n (M) a 
monotonically increasing function of M when w > — 1 and a monotonically decreasing 



approximation scheme remains applicable. 
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Figure 4. Same as figure 3, but with the sound speed fixed at c 2 s = 10~ 6 . 



function of M when w < — 1. The M-dependence is however quite weak for those cases 
with w close to — 1, since dark energy clustering is generally suppressed by a factor 
(1 + w). 

The most interesting case presented here is that for w = —0.8, especially for 
c\ = 10~ 5 and c 2 s = 10" 6 (corresponding Jeans masses: 1.3 x 1O 16 M and 4 x 1O 14 M ). 
Here, we see that at M < Mj, 5 l ™ n is at its lowest value and is essentially independent 
of M, indicating that we are in the fully non-clustering regime. As we move to higher 
values of M, we encounter a transition region where rises with M. Once M 3> Mj, 
however, 8 l ™ n reaches a plateau, where clustering is most efficient and 6^ is again 
independent of M. Interestingly, a similar pattern can also be seen in the w < — 1 cases, 
where the dark energy density contrasts are negative, corresponding to dark energy 
underdensities, which have a negative effect on the clustering of matter. This negative 
effect is strongest for masses larger than the Jeans mass. 

The dependence of 5 l ™ n on the collapse redshift is also quite clear in figure 5. The 
later the collapse, the larger the difference between the values of 5 l ™ n for the clustering 
and the non-clustering solutions. The reason is simply that the contribution of dark 
energy to the total energy budget in the universe increases with time. This means the 
effect of dark energy clustering also becomes more important at lower redshifts. 

Finally, we caution the reader again that our quasi- nonlinear approximation breaks 
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Figure 5. Linear threshold density S 1 ™^ as a function of the halo mass M at various 
collapse redshifts. Top left: Collapse redshift of z = 0. Solid lines denote models with 
w = —0.8 and, from top to bottom, c 2 = 10~ 6 , 10~ 5 , 10 -4 , 10~ 2 . Dashed lines denote 
models with w = —0.9 and c 2 = 10~ 6 , 10~ 2 . For those cases in which 8 1 q > 1 during 
the collapse, we indicate the results with dotted lines. The blue dot-dash and the red 
dot-dot-dot-dash lines represent the reference EdS and ACDM models respectively. 
Top right: Same as top left, but with the solid lines denoting models with w = —1.2 
and, from top to bottom, c 2 = 10 -2 , 10~ 4 , 10~ 5 , 10~ 6 . Dashed lines represent models 
with w = —1.1 and c 2 = 10~ 2 ,10~ 6 . Middle: Same as top panel, but for a collapse 
redshift of z = 1. Bottom: Same as top panel, but for a collapse redshift of z = 2. 



down if the dark energy density contrast inside the top hat becomes too large, especially 
when the halo mass approaches or exceeds the Jeans mass associated with a chosen sound 
speed. As a rule of thumb, we take the condition of breakdown to be 5^ > 1 at any time 
during the collapse process. To alert the reader to those cases where the breakdown 
condition is met, we indicate the resulting linear threshold densities <5j.™ u in figure 5 with 
dotted lines; in these cases, the exact values of S 1 ™^ are unreliable. Clearly, for equation 
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of state parameters that deviate significantly from w — — 1 (e.g., w = —0.8, —1.2), our 
approximate approach breaks down already at M < Mj. For w = —0.9,-1.1, the 
approach appears to remain valid for a larger mass range. 

The breakdown of our approximation in the clustering limit also explains why we 
do not recover the exact results of reference [31] for c\ = and w = —0.7,-1.3. In 
fact, in the clustering limit, our approximation appears to overestimate the effect of 
dark energy clustering on the linear threshold density; had we included all nonlinear 
effects, nonlinear dark energy clustering would feed back on the matter clustering more 
effectively, thereby leading to an earlier collapse for w > — 1 (and a later collapse for 
w < —1). An earlier collapse means that the linearly evolved matter density contrast 
would reach a lower value at the time of collapse, so that the real 8 l ™ n in the clustering 
limit would be lower than our estimate in figure 5. 

4-3. Virialisation 

In reality the collapse of an overdense region will never take place in the way described 
above, since density fluctuations inside the region will moderate the infall, and the 
system reaches virial equilibrium before the matter density can ever become infinite. For 
a single component system such as the case of a dark matter-dominated EdS universe, 
the process of virialisation and the radius at which virial equilibrium is attained can be 
obtained directly from the spherical collapse physics by assuming energy conservation 
between the time of turnaround — defined as the moment at which the top hat radius 
begins to shrink — and the time at which virial equilibrium is established. The result 
turns out to be rather simple, 



i.e., virialisation is complete by the time the top hat radius decreases to half its value 
at turnaround. 

For more complicated systems, such as the two fluid system considered in this paper, 
the conditions of energy conservation need to be modified. The time of virialisation can 
still be taken to be the instant at which the virial theorem is satisfied. However, the 
problem is complicated by the fact that we do not know how or if the dark energy fluid 
takes part in the virialisation process (see, e.g., [47,48]). Therefore, for simplicity, we 
adopt equation (4.6) as the virialisation condition for this paper, and define the virial 
overdensity as 



where t v - 1v is the time at which R = R vir . We expect the quantitative results to be 
somewhat sensitive to our choice of the virialisation condition, but the qualitative 
features should be unaffected. 

Figure 6 shows A vir as a function of the halo mass M for various combinations of 
w and c 2 s . The initial top hat matter density contrast 5^ has been fixed so that all 



vir -^turnaround 




(4.6) 



A 




(4.7) 
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Figure 6. Virial overdensity A v ir as a function of the halo mass M for different dark 
energy equation of state parameters and sound speeds. In all cases 5^ i has been chosen 
such that the collapse occurs at z = 0. Top left: w = —0.8 and, from top to bottom, 
c 2 = l(T 2 ,lCr 4 ,l(r 5 ,l(r 6 . Top right w = -1.2 and c 2 = 1(T 6 , 10~ 5 , 10~ 4 , 10~ 2 . 
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halos collapse at z — 0. As a reference point, for halos that collapse today, A v j r = 147 
and A vir = 252 for an EdS and a ACDM universe respectively (we do not plot these in 
figure 6 because they fall out of the plotting range). 

Similar to the linear threshold density <5'™ u , A vir is first and foremost dependent on 
the choice of w. Introducing a finite dark energy sound speed into the picture induces 
for A v i r a dependence on the halo mass M. However, while ^oii increases with M 
for cosmologies with w > — 1, A vir decreases with it. The opposite trend is seen for 
cosmologies with w < — 1. As with 5'™ u , again we see an asymptotic non-clustering 
value for A vir at M < Mj, a transition region at M ~ Mj where A vir varies strongly 
with M, and a second asymptotic region in the M 3> Mj clustering limit. 

The results in figure 6 are for a collapse redshift of z — 0. For halos that collapse 
earlier, the dependence of A vir on the halo mass is qualitatively similar to the z = case, 
but the difference in A V i r between the clustering and the non-clustering limits is smaller. 
This trend is reminiscent of the results in figure 5 for the linear threshold density 
for different collapse redshifts. 

Finally, note that some authors define the virial overdensity as the top hat density 
at the time of virialisation r vir , but normalised to the background density evaluated at 
the collapse time r co ii. This means that instead of, e.g., A v ; r = 147 for the EdS model 
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according to our definition (4.7), one finds a higher value of 179 simply because between 
T vir and r co ii the background density has become smaller due to the Hubble expansion. 
Since the dark energy sound speed does not play a role in the background expansion, 
the effect of these differing definitions is only to induce a shift in the normalisation of 
A vir for a given set of w and VLq. The M-dependence of A vir is unaffected. 

5. Discussions and conclusions 

While a dark energy fluid with a negative equation of state parameter appears to describe 
the apparent accelerated expansion of our universe with reasonable success, the precise 
nature and properties of this dark energy remain undetermined. In this paper, we 
have addressed some aspects of the dark energy's role in cosmic structure formation. 
Specifically, we have considered a generic dark energy fluid parameterised by a constant 
equation of state parameter w and sound speed c s , and determined their impact on the 
formation of gravitationally bound objects. 

Our main tool is the spherical collapse model, incorporating a nonrelativistic dark 
matter component and a generic dark energy fluid described above. Such a model has 
been investigated by other authors previously in the limit where the dark energy is (i) 
non-clustering, i.e., c s — > oo, or (ii) comoving with the dark matter, i.e., c s — > 0. In this 
work, we have generalised the spherical collapse model to describe those intermediate 
cases characterised by an arbitrary c s . Although we have focussed on scenarios with 
constant w and c s , our approach is equally applicable to dark energy fluids described by 
time-dependent parameters. Along the way, we have also provided a detailed description 
of dark energy evolution in the linear regime — again for arbitrary w and c s , and identified 
some salient features of dark energy clustering. 

The dark energy component changes the evolution of the spherical collapse through 
its effect on the overall expansion of the universe as well as through its own clustering 
abilities. In the non-clustering limit, only the equation of state parameter w plays a role 
in the dynamics of the spherical collapse by altering the rate of the Hubble expansion. 
When the dark energy is able to cluster, however, then the dark energy density contrast 
also sources the evolution of the dark matter density perturbations. The amount of 
dark energy clustering is determined by both w and c s . 

In addition, since the introduction of a dark energy sound speed c s necessarily brings 
into the picture a "Jeans scale" and hence a "Jeans mass" Mj, we find that the same 
sound speed can influence the spherical collapse dynamics in different ways depending 
on the mass of the collapsed object or the "halo mass" M. This mass dependence 
is especially manifest when we compute the virial overdensity A v i r and the linearly 
extrapolated threshold density 5^ u . In both cases, we find two asymptotic regions 
corresponding to the non-clustering M <C Mj limit and the clustering M 3> Mj limit, 
where the value of A vir or 5 l ™ n is practically constant with respect to M. In between 
these limits is a transition region, in which A vir (6^) declines (grows) strongly with M. 
Observing this transition region will be tell-tale sign that dark energy is dynamic, and 
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a great leap towards pinning down its clustering properties. 

One possible way to discern this halo mass- dependent dark energy clustering is to 
measure the cluster mass function. Already it has been shown in, e.g., reference [31] 
that the difference between a clustering [c s = 0) and a homogeneous (c s — > oo) dark 
energy fluid can be an order unity effect on the expected number of clusters at the high 
mass tail of the mass function. It remains to be seen how exactly an arbitrary dark 
energy sound speed would alter this conclusion, but some additional scale-dependent 
effects are almost guaranteed to be present. 

Finally, let us remind the reader again that in our spherical collapse analysis we 
have simplified the equations of motion for the dark energy component so that only 
terms linear in the dark energy density contrast have been retained. This approach is 
strictly not valid for those cases in which the dark energy density contrast exceeds unity 
during the collapse process, and is especially prone to breakdown for cosmologies with 
an equation of state parameter that deviates significantly from w = —1. In these cases, 
a fully nonlinear analysis is required in order to compute accurately such quantities 
as A vir and 8 l ™ n for comparison with observations. However, the qualitative features 
of dark energy clustering should not be affected by our simplified approach. A fully 
nonlinear analysis would require that we solve the equations of motion for the dark 
energy component — two (1 + l)-dimensional partial differential equations — using either 
a grid-based finite difference scheme or a Lagrangian method akin to a smoothed particle 
hydrodynamics simulation. Investigation is already underway and we hope to report the 
results in a future publication. 
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Appendix A. From general relativity to the pseudo-Newtonian approach 

We demonstrate in this section that, at linear order, the pseudo-Newtonian approach 
adopted in this work is indeed consistent with the Newtonian limit of a general 
relativistic formulation. We work in the conformal Newtonian gauge, in which the 
perturbed line element is given by [49] 

ds 2 = a 2 (r){-[l + 2ip(r, x)]dr 2 + [1 - 20(r, x)}(dx 2 + dy 2 + dz 2 )}, (A.l) 

for a flat background spatial geometry. Assuming zero anisotropic stress so that 
i/j — (j), the equations of motion in Fourier space for the dark matter and dark energy 
components [37-39] are, respectively, 

S m + 9 m - 3<p = 0, 

L + W m - k 2 <p = 0, (A.2) 
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1 + w 

where D = l + 9(c 2 s — w)l-i 2 /k 2 , and D = ld 1 {c 2 — w)'H'H/k 2 . Rearranging equations (A. 2) 
and (A.3) respectively into second order differential equations, we find 

L + U8 m = 3H<j) - k 2 <f) + 30, 



(A.4) 



and 



S Q + [(1 - 3w)H - D/D]5 Q 

+ {Dc 2 k 2 + 3(c 2 -w)[H+(l- 3c 2 s )n 2 - (D/D)n]}5 Q 

= (1 + w){3[(l - 3c 2 s )H - D/D]<j) - Dk 2 <f) + 30}, (A.5) 

assuming time- independent w and c 2 . 

In order to take the subhorizon (i.e., k 3> "H) limit of equation (A.5), we note that 
"H ~ 0("H 2 ). Thus, we can replace all occurrences of D with D = 1 and D/D with 
D/D = 0. One last step concerns the coefficient of the Sq term on the LHS: we group 
all contributions proportional to c 2 together to get 

.w 2 (n 
\w + 

and set the 0("H 2 /k 2 ) terms in [. . .] to zero. Thus, we find 

6 Q + (1 - 3w)W Q + [c 2 k 2 - 3w(H + n 2 )]5 Q 

= (l+w) [3(1 - 3c 2 )^0 - k 2 (f) + 30] (A.7) 

for the subhorizon limit of equation (A.5). 

It remains to specify the metric perturbation <fi in terms of the density and velocity 
perturbations via the Einstein equation. Using the expressions given in, e.g, [49], we 
obtain 



2i.2 



cik 



1 + 3 



3(6 



w) 



3w(U 2 + H) \ S Q , 



(A.6) 
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For the problem at hand, the summation is performed over a = m, Q, and the notation 
should be understood to mean Q a = Q a (r), w rn = 0, wq = w, c 2 m = 0, and Cq = c 2 . 
Combining these expressions to form the RHS of equation (A.4), we find the leading 
order contribution in the subhorizon limit to be 

3-H0 - k 2 <p + 30 ~ -H 2 EO*(l + Zc 2 a )5 a . (A.9) 
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Similarly, we find for the RHS of equation (A. 5), 

3(1 - 3c 2 s )U<P - k 2 <p + 30 ~ \u 2 ^fi a (l + 3c 2 a )6 a , (A.10) 

again to leading order mTi/k. Note that RHS of equations (A. 9) and (A. 10) are identical 
as a result of a % 2 //c 2 -suppressed l-Lcf) term relative to the k 2 <ft and <j) terms, as can be 
seen in equation (A. 8). We have not assumed c 2 <C 1 to arrive at this result. 
Thus, equations (A. 4) and (A. 5) now become, respectively, 

'L + nL = ^n 2 [n m (r)5 m + n (r)(i + 3c 2 s )5 Q ], 

S Q + {1- 3w)U5 q + [c 2 s k 2 - 3w(n + n 2 )}5 Q 

= (1 + w)^n 2 [n m (r)S m + n Q (r)(l + 3c 2 )5 Q }. (All) 

These equations are consistent with the outcome of the pseudo-Newtonian approach 
(equations (3.1) to (3.3)). 

Appendix B. Validity of the H 2 S > k 2 > H 2 regime 

A number of approximations have been made in order to arrive at the final equation 
of motion (A. 11) for the dark energy component in the subhorizon limit. We have 
demanded that (1 — 3u>)"H ^> D/D in equation (A. 5), which limits the use of the 
approximate equation (A. 11) to H 2 /k 2 0.074 — > 0.15, depending on the sound 
speed assumed. We have also assumed D ?a 1, equivalent to imposing a limit of 
l-L 2 /k 2 0.055 — > 0.1 on the validity of our approximate equations. In practice, we 
might want to choose a benchmark limit of 

max (^) = 0.01, (B.l) 

in order to keep all superhorizon contributions at a truly subdominant level. 

Similarly, if we wish to take the super-sound-horizon limit (i.e., the clustering limit, 
k <^ % S i where % s = l-L/c s ), then from equation (3.2) or (A. 11), we see that 

maX (^) = °- 1 (R2) 
would work as a good benchmark limit to ensure the subdominance of sub-sound- 
horizon contributions. The two benchmark limits (B.l) and (B.2) together define the 
T-L 2 S ^> k 2 ^> H 2 regime (or the "clustering regime" in this work). 

From here, it is easy to see that the two limits (B.l) and (B.2) combine to set a 
constraint on the dark energy sound speed of 

c\ < lO" 3 (B.3) 

in the H 2 S> k 2 3> H 2 regime. If we were to exceed this constraint, we also run the risk of 
increasing the superhorizon ("H 2 /fc 2 ) and/or sub-sound-horizon (k 2 c 2 /'H 2 ) contributions 
to any approximate analytic solution to beyond the "negligible" level, in which case the 
1-L 2 ^>k 2 ^> Ti 2 regime becomes ill-defined. 
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